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We describe a number conserving approach to the dynamics of Bose-Einstein condensed dliute atomic gases. 
This builds upon the works of Gardiner [C. W. Gardiner, Phys. Rev. A 56, 1414 (1997)], and Castin and Dum 
[Y. Castin and R. Dum, Phys. Rev. A 57, 3008 (1998)]. We consider what is effectively an expansion in powers 
of the ratio of non-condensate to condensate particle numbers, rather than inverse powers of the total number 
of particles. This requires the number of condensate particles to be a majority, but not necessarily almost 
equal to the total number of particles in the system. We argue that a second-order treatment of the relevant 
dynamical equations of motion is the minimum order necessary to provide consistent coupled condensate and 
non-condensate number dynamics for a finite total number of particles, and show that such a second-order 
treatment is provided by a suitably generalized Gross-Pitaevskii equation, coupled to the Castin-Dum number- 
conserving formulation of the Bogoliubov-de Gennes equations. The necessary equations of motion can be 
generated from an approximate third-order Hamiltonian, which effectively reduces to second order in the steady 
state. Such a treatment as described here is suitable for dynamics occurring at finite temperature, where there 
is a significant non-condensate fraction from the outset, or dynamics leading to dynamical instabilities, where 
depletion of the condensate can also lead to a significant non-condensate fraction, even if the non-condensate 
fraction is initially negligible. 

PACS numbers: 03.75.Nt, 67.40.Db, 05.30.Jp 



I. INTRODUCTION 



By definition, a dilute atomic gas that has undergone Bose- 
Einstein condensation |lj, |2j, |3[] has a large number of compo- 
nent particles occupying the same mode [H H H H 0, H Ml ■ 
Effects associated with such a macroscopic occupation were 
first observed in superfluid helium and in superconducting 
metals ]9|]. The importance of interactions in such com- 
paratively dense condensed-matter systems means that the 
condensate fraction, although important, is substantially less 
than the non-condensate fraction. In systems composed of 
laser and magnetically cooled and trapped dilute atomic gases 
lfl0l[TTl[l2ll the situation is often very different; the atomic gas 
can be sufficiently cold and dilute for the condensate fraction 
to be a large proportion of the total number of atoms. It is for 
this reason that the Gross-Pitaevskii equation mHEHO, 
originally conceived to develop a qualitative understanding of 
processes in superfluid helium, has achieved the status of a 
quantitatively useful description of degenerate dilute gases of 
bosonic atoms. 

The Gross-Pitaevskii equation is essentially a classical field 
approximation to an underlying quantum field. Notwithstand- 
ing its broad utility, there are many situations where a more 
accurate description is required. Superfluid to Mott-insulator 
phase-transitions in optical lattices 017I1 . and dimer formation 
via controlled manipulation of magnetic fields (in order to ex- 



ploit Feshbach resonances) lfl8ll are topical examples of such 
processes. Effectively the strength of the inter-atomic interac- 
tions becomes significant to the extent that higher-order atom- 
atom correlations must be more carefully accounted for, and 
for which the standard Gross-Pitaevskii equation is inadequate 
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Even apart from such extreme situations, if the non- 
condensate fraction becomes significant, a description going 
beyond the Gross-Pitaevskii equation must be called upon. 
Two important situations where this may occur are: dynamics 
occurring at a (significant) finite temperature 025I1 . of interest 
due to the unique possibility offered by dilute Bose-Einstein 
condensate experiments for quantitative tests of thermal field 
theories; and dynamics leading to dynamical instabilities in, 
and hence depletion of an initially low temperature conden- 
sate ll26i l27l l28ll . such as may well occur in experiments 
ll29l, l30l l3lll studying chaotic an d q uantum chaotic dynam- 
ics in Bose-Einstein condensates ll32l l33i l34l l35l l3rjl 13711 . The 
desire to provide a relatively simple, consistent description of 
condensate and non-condensate dynamics motivates the work 
presented here, and a form of the approach we present was a 
key part in work carried out l38l HsJ H^], to good agreement 
with experiment B25I) . in order to describe excitations at finite 
temperature of a dilute Bose-condensed gas. 

The first recourse when wishing to go be yon d the Gross- 
Pitaevskii emationilJlJMMSl 
IS gl |M IB |H ElMMMMis frequently the Bo- 
goliubov, or Bogoliubov-de Gennes equation s ll4l[ l42i l43ll . 
or their number-conserving variants 11441 l45l 14a, 14711 . Par- 
ticularly motivated by the desire to explain the properties 
of Bose-condensed gases at finite temperature, a number of 
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extensions have been proposed. These include generaliza- 
tionsj5ipl Hi E>H lH IH B of linear response the- 
ory stochastic nterpretations of the Gross-Pitaevskii 
equation 1671 roa, ro9l p70l 17111. Hartree-Fock-Bogoliubov ap- 
"IJ^^lZlJlRJl^JMLa varietyof 



proaches 

kinetic theories [81, 82, sllMllSjMM 



and a cumulant-based formalism 112 lL l93l 

The description presented here is within a number- 
conserving formalism, and builds on the works of Gardiner 
Hil], and Castin and Dum l45ll . which are essentially equiv- 
alent to each other. Symmetry-breaking formulations, which 
automatically violate particle number conservation, have met 
with considerable success in describing the observed proper- 
ties of Bose-Einstein condensed dilute atomic gases. Here 
symmetry-breaking is defined as the breaking of the U(l) 
global phase symmetry whose Noether charge is the total par- 
ticle number. Technically, however, symmetry-breaking for- 
mulations require a coherent superposition of different num- 
bers of particles. One could argue that the actual particle 
number is only known statistically in any real experiments, 
and should be considered an ensemble average from multi- 
ple realizations of the same experiment. Even given this, it 
is difficult to see how shot-to-shot number-coherences could 
be built up. It is therefore important to understand any dif- 
ferences which might appear between number-conserving and 
symmetry-breaking formulations. 

A significant difference that does occur, in the number- 
conserving formulation presented in this paper, is the presence 
of nonlocal terms in the equations of motion for the conden- 
sate and non-condensate fractions. These arise from defining 
the two fractions in such a way that they must be mutually or- 
thogonal. This orthogonality is not in general fulfilled if one 
defines the condensate fraction as being the expectation value 
of a field operator, as occurs in symmetry-breaking formula- 
tions. The nonlocal terms can have a significant effect; correct 
inclusion of these terms has been observed to be crucial in ob- 
taining good agreement between application of the number- 
conserving formulation presented in this paper l38l l39l l4(ill . 
and the experimental study of collective excitations in a finite- 
temperature Bose-Einstein condensate l25ll . 

We consider what is effectively an expansion in powers 
of the ratio of non-condensate to condensate particle num- 
bers, rather than inverse powers of the total number of par- 
ticles [44, 45]. This means that there should be signifi- 
cantly more condensate than non-condensate particles, but 
that the condensate fraction is not necessarily assumed to be 
nearly encompassing the entire many-body system. We ar- 
gue that a second-order treatment (in the dynamical equa- 
tions of motion) is the minimum order necessary to pro- 
vide consistent condensate and non-condensate number dy- 
namics, with exchange of particles between the fractions, 
for a finite total number of particles. We show that such 
a second-order treatment is provided by a suitably general- 
ized Gross-Pitaevskii equation, coupled to the Castin-Dum 
number-conserving formulation of the Bogoliubov-de Gennes 
equations (these are modified only by the presence of projec- 
tors necessary to maintain orthogonality between the conden- 
sate and non-condensate components) [45]. The necessary 



equations of motion can be generated from an approximate 
third-order Hamiltonian, which effectively reduces to second 
order in the steady state. 

This paper is organized as follows: Section HTl formally de- 
scribes the many-Boson system under consideration, and de- 
termines a suitable fluctuation operator on which to base the 
expansion; Section|III]constructs an appropriate cubic approx- 
imate Hamiltonian used to generate the desired equations of 
motion, and justifies the approximations made; Section [IV] 
elucidates the equations of motion detailing both condensate 
and non-condensate dynamics, systematically to zeroth, first, 
and second order in the fluctuation operators; Section [V] dis- 
cusses some considerations when the system is assumed to be 
in an equilibrium state; and Section[VT]consists of the conclu- 
sions. There then follow two technical appendices elaborating 
on points made in the main text, which are included for com- 
pleteness. 



II. SYSTEM PROPERTIES 
A. Model Hamiltonian 

The starting point of the theory is the binary interaction 
Hamiltonian for a system of bosons, 



H(t) 



I 



<W (r)// sp (r, f)*P(r) 



"if 



(1) 



c/r»P T (r)T' 1 (r) > P(r)T'(r). 



The field operators obey the standard bosonic commutation 
relations [*P(r), ^(r')] = S(r - r'), and 

h 2 



tfsn(r) 



2m 



V z + V(r) 



(2) 



(where m is the atomic mass) is the single-particle Hamilto- 
nian, containing the kinetic energy and any external potential 
energy terms. 

For simplicity we have assumed the binary interaction to 
be characterized by an energy-independent contact potential, 
Vbin(r - r') = UoS(r - r'), where Uq = 4-7th 2 a s /m and a s is 
the s-wave scattering length. This is a standard approxima- 
tion for three-dimensional, cold dilute Bose gases. As is well 
known, however, it leads to ultra-violet divergences, which 
are removed by renormalizing various quantities appearing in 
the subsequent development of the theory. This procedure is 
well understood and has been discussed by a number of au- 
thors (see, for example, Refs. lEH Izl, I^IL H^, H3, HM] ) • It 
can be rigorously justified, and we give a brief outline of the 
relevant arguments in AppendixlAl 



B. Condensate and fluctuation terms 

1. Single-body density matrix 

In the same manner as Castin and Dum l45ll . we follow 
Penrose and Onsager 19911 in defining the condensate wave- 
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function. This is in terms of the single-body density matrix 
p(r,r',0 = <*VyF(r)>. (3 ) 

The single-body density matrix is Hermitian, i.e., p(r, r', if = 
p(r',r, f), and can be decomposed into a complete set of 
eigenfunctions with real eigenvalues. As p(r, r', t) may be 
time-dependent, these are the instantaneous eigenfunctions 
and eigenvalues, defined by the diagonal representation of the 
single-body density matrix at a specific time t . 

We assume that there is one distinct eigenfunction 0(r, f), 
defined with unit norm, which has a corresponding eigenvalue 
N c (t) significantly larger than all the other eigenvalues. Thus 

j dr'p(r, r', f)cf>(r', t) = N c (t)(f>(r, t) (4) 

(time arguments are used to indicate a possible explicit time 
dependence). We are free to partition the field operator into 
condensate and non-condensate components: 

*(r) = a c (t)<f>(r, t) + 5%{r, t), (5) 

where a c {t) annihilates a particle in mode </>(r, f), and <J*F(r, t) 
is that part of the field operator 4*(r) orthogonal to (f>(r, t). We 
refer to 0(r, f) as the condensate mode. 

2. Comparison with symmetry-breaking 

The conventional symmetry-breaking methodology to par- 
tition the field operator into condensate and non-condensate 
components assumes a finite expectation value *F(r, f) = 
(4*(r)) for the field operator, and additionally considers a fluc- 
tuation term <s4* s b(r, f) around this mean value JH S, 0, S- 
Hence, one states: 

*(r) = 'P(r,t)+a*.b(r,0. (6) 

By definition, a fluctuation around a mean or expectation 
value must itself have a mean or expectation value equal to 
zero, i.e., (^^(r, f)) = 0. It also follows that these fluctu- 
ation operators have exactly bosonic commutation relations, 
i.e., 

[d* lb (r, f), SV A <J, t)] = [%), 4V)] = S(r - r'). (7) 

Comparing the partitioned terms arising from a symmetry- 
breaking formulation with those arising from consideration of 
the single-body density matrix, we get 

«P(r, t) = <fi c (f»0(r, t) + <J*(r, f)>, (8) 
S# sb (r, f) =[a c (0 - (a c (t))]<f>(r, f) 

To enable a direct comparison with symmetry-breaking, we 
have not assumed particle-number conservation, and hence 
neither (a c {t)) nor (5*T(r, f)) are automatically assumed = 0. It 
is apparent that exactly which parts of the many-body system 
are designated as condensate and non-condensate, depends 



on the formulation employed. One consequence is that, un- 
like </>(r, t) and <54*(r, f), the terms ^(r, t) and <j4* s b(r, t) aris- 
ing from a symmetry -breaking formulation are not, in general, 
spatially orthogonal. As can be seen from Eq. ([8]l and Eq. (O, 
in general they both have contributions from both 0(r, t) and 
<5*F(r, t). 

From the Penrose-Onsager Ir99ll formulation for partitioning 
the field operator used in this paper, it follows that those parts 
of the many-body system designated as being condensate and 
non-condensate must always be orthogonal. It is this require- 
ment that leads to higher-order (i.e., beyond Gross-Pitaevskii) 
equations of motion necessarily including nonlocal terms, not 
arising in symmetry-breaking formulations. It should be em- 
phasized that this difference does not in itself imply that one 
formulation is more correct than the other, although partition- 
ing the field operator into manifestly orthogonal components 
can be a considerable convenience. 



3. Commutation relations 

Formally, the condensate mode creation operator a\(t) and 
the non-condensate field operator ^(r, f) can be defined with 
respect to the bosonic field operator 4*(r): 

aj(f) = j dr¥(r)(f>(r,t), (10) 

where the projector Q(r, r', f) is defined to be 

Q(r, r', f) = S(r - r') - 0(r, t)f(r\ t). (12) 

It follows that the only non-zero commutation relations in- 
volving a c ( f), S^ir, f), and their Hermitian conjugates are: 

[a c (t),tl{t)-\ = \, (13) 
[^(r, f), SV\r', f)] = Q(r, r', t). (14) 

C. Fluctuation expansion 

1. Overview 

It is a common procedure to consider a mean value of some 
observable quantity, and to account for the observed variance 
in the values obtained from multiple runs of some real or 
hypo thetical measurement procedure with a fluctuation term 
lllOOll . The observable quantity is then described in terms of 
its mean or expectation value, plus a fluctuation term with (by 
definition) zero expectation value. If the observable is a sim- 
ple one-dimensional quantity, such as the vertical position of a 
classical point-mass, then the statistics of the measured values 
can be represented by a simple one-dimensional distribution. 

If the fluctuations are vanishingly small, then the distribu- 
tion of a continuous variable collapses to a 5-function, but 
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otherwise the distribution will have a width, which one usu- 
ally describes in terms of the second-order cumulant (the vari- 
ance). The scale of the fluctuations is therefore set by the 
square root of the variance (the standard deviation). Con- 
versely, finite fluctuations imply a finite width, and hence a fi- 
nite variance. In the case of a Gaussian distribution, all higher- 
order moments and cumulants can be described in terms of the 
distribution's mean and variance. 

The symmetry-breaking methodology used to partition the 
field operator, as described in Section HlB 21 is an example of 
describing an operator in terms of its expectation value plus 
a fluctuation term. Analo gues to the variance are the second- 
order cumulants ll93l flOOh : 

(❖Vmr')) - nr, fW, f) = <<^ b (r, f)^ s b(r', f)>, (15) 
»')) - »F(r, W, t) = <5* sb (r, 0(S*sb(r\ f)>, (16) 

together with their Hermitian conjugates. As 
( j dr6^\(r, f)<J*i'sb(r, 0) is the number of non-condensate 
particles, we generally consider 5x s i,(r, t) to scale as the 
square root of the number of non-condensate particles. Note, 
however, that even if all second- (and higher-) order normally 
ordered cumulants equal zero, there still remain purely 
quantal fluctuation terms, due to the noncommutativity of the 
field operators. 

If the physically relevant matrix elements of the fluctuation 
terms can be considered to be in some sense small (for exam- 
ple with respect to the square root of the number of condensate 
particles, as considered in this paper), then one can attempt an 
expansion of the Hamiltonian [Eq. (Q])] and equations of mo- 
tion in terms of products of the fluctuation terms. Terms of 
beyond a specified order can then be discarded as negligible, 
up to the order of interest S S 0, IE E2] - 

In this paper, we wish to carry out a comparable program, 
but within a number-conserving formalism li44l l45l l57ll . As 
the expectation values of all particle creation and annihila- 
tion operators are (trivially) equal to zero, we cannot define 
a "small" fluctuation operator by taking a field operator and 
subtracting its (zero) expectation value. An appropriate fluc- 
tuation operator that does not trivially have zero expectation 
value is required. The validity of such expansions is in general 
reliant upon (at least in the homogeneous limit) (Na\) 1 ^ 2 <K 1 
if T = 0, and (k b T /N c U )(N c al) 1/2 « 1 if (k b T/N c U ) » 1, 
where N is the total particle number, T is the temperature and 
ks is Boltzmann's constant lf57ll . 



2. Number-conserving fluctuation operators 

Substituting Eq. © into Eq. © reveals that N c (t) = (N c (t)), 
where N c (t) = ci\(t)a c {t). The eigenvalue N c (t) is thus the 
mean number of particles in the condensate mode. Further- 
more, 



(al(f)SV(r, f)> = 0, 



(17) 



i.e., there are no simple coherences between the condensate 
and non-condensate components. 



The product d£(t)5¥(r,t) has much to recommend it as a 
fluctuation operator suited to a number-conserving formal- 
ism. Unlike ^(r, f), its mean value is not trivially equal to 
zero in a number-conserving approach. The desired number- 
conserving fluctuation operator should be of the same magni- 
tude as 6^(r, t), however. 

Castin and Dum l45ll . and Gardiner l44tl were thus moti- 
vated to define 



A(r,f) 



1 



N 



aj(o^P(r, o, 



(18) 



where A^ = f dr^(ry9f(r) is the total particle number op- 
erator. Introducing the notation N,(t) = N - N c (t) for the 
non-condensate fraction, and noting that ct c (t) should scale 
as VNc(0 and <54*(r, t) as \W,(f), it then follows that A(r, t) 
scales as 



N c {t)N,(t) 



N 



N,(t) 



N,(t) 



N 



(19) 



for N » N,(t). Hence, under the assumption that N c (t) a N, 
this operator scales satisfactorily. 

When considering an assembly of exactly N atoms, Eq. ( TPTb 
directly implies that (A(r, f)) = 0, and so A(r, t) can be con- 
sidered a simple fluctuation operator. To linear order in A(r, t) 
JH, 

[A(r, t), A + (r', f)] ~ [W(r, 0, ^ f (r', t)] = Q(r, r', f), (20) 
and 



drA\r, t)A(r, t) 



I 



drtf¥\r, t)8¥(r, t) = N - N c (t). 

(21) 

Again, when considering an assembly of exactly atoms, 
there can be no fluctuations in the total number operator, and 
so N - N c (t) can be identified with N - N c (t). 

We wish to avoid making the assumption that Af c (f) ~ N, 
i.e., that almost all bosons are in the condensate mode, and so 
consider a scaling proportionate to the number of condensate 
atoms rather than the total number of atoms J^l4^ll47[l57ll . 

A possible alternative to A(r, t) [Eq. ( fT8l )l is 



A c (r,f) = 



1 



zal(f)6W(r, t), 



(22) 



from which the exact identities 

[A c (r, f), Aj(r', f)] = [d*(r, f), SWtf, t)] = Q(r, r', t), (23) 
and 

J drkl(r, t)A c (r, t) = j dr6^(r, f)S¥(r, t) = N - N c (t). 

(24) 

follow. This strong correspondence between normal pairs of 
A c (r, t) and 6^(r, t) operators appears very attractive. How- 
ever, the expectation value (A c (r, f)) is not guaranteed to be 
identically equal to zero. Consequently, the operator A f (r, t) 



5 



cannot be treated as a simple operator- valued fluctuation. This 
complicates the development of a consistent expansion in 
terms of products of A c (r, t) for the determination of improved 
equations of motion. 

In particular, the methodology introduced by Castin and 
Dum [45], and employed in slightly modified form in Sec- 
tion [TV] of this paper, relies explicitly on {dA(r, f)/dt) = 
d(A(r, f))/dt = for the development of the equations of mo- 
tion. As (A c (r, t)) is not exactly equal to zero, an equivalent 
expression is not guaranteed to hold for A c (r, t). One could 
still exploit the idea that (A c (r, f)) should be approximately 
equal to zero, but this would in principle require a careful ac- 
counting of higher-order corrections to such an approxima- 
tion, for example involving Taylor expanding [iV c (0] -1 ' 2 m 
powers of [N c (t) - N c (f)]/N c (t). 

Equation ( TTTb tells us that al(i)SY(r, t) and any scalar mul- 
tiple thereof has expectation value exactly equal to zero. In 
order to avoid such complications as described above, we thus 
choose to carry out an expansion in terms of 

A(r,f) = —^=at(f)6^(r,f). (25) 

The normal A(r, t) pair is related to the normal 6^{r, t) pair 
via 

A'Xr', f)A(r, t) = ~^r^ 1 SW</* t)S$(x, ?), (26) 
and the exact commutation relation is given by 

[A(r,f),AV,0] =^e(r,r',f) 

-^tf# t (r , ,f)d*(r,0. 

3. Comparison of the candidate fluctuation operators 

The operators A(r, t) and A c (r, f), in comparison to A(r, t), 
scale more straightforwardly as y/N t (t). In the expansion 
used in this paper, the terms that appear in the Hamilto- 
nian and equations of motion have the form of products of 
A(r, t)l y/N c (t), and/or its Hermitian conjugate [there is an 
overall factor of N c (f) in the Hamiltonian, and an overall factor 
of y/N c (t) in the equations of motion for A(r, t) and A 1 (r, f)]. 
The fluctuations must thus be small compared to the square 
root of the number of condensate particles. In other words, the 
small parameter such a fluctuation expansion is y/N,(t)/N c (t) 
l45ll . At the very least, greater than N/2 particles must there- 
fore be in the condensate mode for the truncation of higher- 
order terms to be justifiable. If UqN is kept constant, then, at 
zero temperature, N, tends to a constant value as N — > 00 (the 
quantum depletion). Hence, 

Va^T V~ k vf (28) 

and one can consider the small parameter in an expansion in 
terms of A(r, t), A f (r, t) to be 1/ V^V jH]. 



In Sections [TTT] and [IV] we will include terms of or- 
der [1/N c (i)] ' 2 whenever we include terms of order 
[N t (t)/N c (t)] k/2 (k an integer). We are obliged to do this if 
we wish to include the limiting case described by Eq. (l28l in 
the formalism, although, as the non-condensate fraction in- 
creases, the contribution of such terms relative to terms scal- 
ing as [N,(t)/N c (t)] k/2 decreases. 

In Section [VBl we will see that a direct consequence of Eq. 
(l27l > being only approximately equal to [6 } ¥(r,f),5 ]< ¥'(r',f)] 
is that a quasiparticle formulation produces quasiparticle 
creation and annihilation operators that are only approxi- 
mately bosonic. This is equally the case for the commutator 
[A(r,0,A t (r',0] [Eq. ®], but not for [A c (r, t), AJ(r', t)] = 
[6^(r,t),6^'(r',t)] [Eq. A problem that therefore 

arises, not present in comparable symmetry-breaking formu- 
lations, is that the desire to have a well-defined fluctuation op- 
erator does not appear to be fully compatible with the desire 
for this fluctuation operator to induce bosonic quasiparticle 
commutation relations. 

One could take the view that, as non-zero corrections to 
(A c (r, f)) would only appear at a higher order than is being 
considered in this paper, one could equivalently consider an 
expansion in terms of A c (r, t) up to the order of current in- 
terest ll39tl . This effectively erases any difference between 
A c (r, t) and A(r, t), however, and it is more straightforward, 
especially when determining truncations of the many-body 
Hamiltonian necessary to generate the equations of motion, 
to consider an expansion in terms of A(r, f) from the outset. 
In particular, doing this avoids the complication of having to 
approximate square-root condensate number operators to an 
appropriate order in the fluctuation operators. This appears to 
cause particular complications in deriving a correct approxi- 
mate third-order Hamiltonian compatible with the equations 
of motion calculated to the corresponding order (although not 
the approach presented in this paper, the equations of motion 
can be deduced directly, without first determining an approx- 
imate Hamiltonian from which they should be derived, as in 
Ref. I45I1 ). This does leave somewhat open the question of 
what the best approach is if one wishes to extend the theory to 
include higher-order terms JH HHH |H . 

We now summarize the relative merits of A(r, ?), A c (r, t), 
and A(r, t). We wish to have a non-trivial fluctuation opera- 
tor that scales like S^(r, t), i.e., as ^N,(t). Both A c (r, t), and 
A(r, f) have this property, whereas this is only true for A(r, t) 
in the limit N t (f)/N -> 0. Of A c (r, t) and A(r, t), A(r, t) is a 
simple fluctuation operator (expectation value exactly equal to 
zero) but corresponds to not-exactly-bosonic quasiparticle op- 
erators, whereas A c (r, t) corresponds to exactly bosonic quasi- 
particle operators, but does not have expectation value exactly 
equal to zero and is therefore not a simple fluctuation operator. 
In this paper we have made the decision that, for the develop- 
ment of the equations of motion, it is simpler to consider an 
expansion in terms of A(r, t). In particular, this avoids the 
difficulty of having to approximate square-root condensate- 
number operators in a consistent way. 
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4. Fluctuation statistics 



In a similar way to how the first manifestation of a fluc- 
tuation about a real mean value is in the variance of its cor- 
responding distribution, for a finite number of particles, the 
presence of fluctuation operators effectively implies non-zero 
values for such pair expectation values as (A'(r',i)A(r,i)). 
This is, of course, with the important caveat that, as A(r, t) 
and A' (r', t) do not commute, there will always be quantum 
fluctuations even if (A* (r', t)A(r, ?)) = 0. With an interacting 
gas, however, (A 1 (r', t)A(r, t)) always has a finite value. 

We thus insist a priori that equations of motion should be 
taken to quadratic order in products of the fluctuation oper- 
ators A(r, t) and A^r, t). This paper lays out a quite gen- 
eral framework for the development of higher-order theories 
involving both the condensate and non-condensate fractions. 
However, we also have a very specific goal, which is to deter- 
mine a minimally complicated set of equations for the self- 
consistent treatment of condensate and non-condensate dy- 
namics. Bearing this in mind, a straightforward simplifica- 
tion is to enforce that all possible expectation values of prod- 
ucts of the fluctuation operators A(r, t) and A 1 (r, f) are either 
zero (for odd products of fluctuation operators), or expressible 
in products of pair expectation values. This is essentially a 
Gaussian approximation, i.e., one assumes that all cumulants, 
or connected correlation functions, of order greater than two 
can be considered negligible B93I1 . One can in principle extend 
this approximation to also account for the existence of higher- 
order cumulants (and hence higher-order correlations). The 
quantitative validity of such an approximation will depend on 
the specific dynamics under consideration. 

This in turn implies that the the many-body Hamiltonian 
[Eq. ([T])] should be approximated to cubic order in the fluctu- 
ation operators [Eq. (1251) 1, This is the minimum order neces- 
sary to produce equations of motion to quadratic order, and it 
is not our intention in this paper to account for any terms of 
any higher order. 



of A(r). Defining U = U[)N C , the result of carrying this out is: 



n c r 



dr 



0*(r)ff sp (r)0(r) + ^_i2||0( r )| 4 



6*(rMr)| 



Nc—l? 



Va£ J dr[f(r)H s? (r)A(r) + H.c] 

I 



-A(r) + H.c. 



+ rfrA'(r) 



N c N c - 1 - , 

/#sp(r) + — 2t/|0(r)| 2 

N c N c 



A(r) (29) 



J dr [<f(r) 2 A(r) 2 + H.c] 



0*(r)A 1 "(r)^A(r) 2 + H.c. 

N c 



N c J 



drA'Ury 



Ni 



N C (N C - 1) 



A(r) 2 , 



where the terms are arranged in ascending order of products 
of the fluctuation operators A(r) and A' (r). 

Equation ( 1291 is an exact reformulation of Eq. (Q]); note, 
however, that A(r) cannot be straightforwardly expanded in 
terms of exactly bosonic quasiparticle operators (see Section 
IV Bl i. and formulating the many-body Hamiltonian in terms of 
bosonic quasiparticle operators can be of great utility in deter- 
mining, for example, energy spectra to high order in a consis- 
tent fashion [57]. It is relatively straightforward to determine 
an equivalent formulation to Eq. (|29l in terms of A e (r) [Eq. 
(l22lil. although this introduces square-root number-operator 
terms yN~c, which can be awkward to deal with. 

As, in the steady state, the highest-order Hamiltonian 
considered in this paper is effectively only second-order in 
the number-conserving fluctuation operators A(r), A T (r) (see 
Section HV D 4K in the present context such considerations can 
largely be avoided. 



B. Reduction to a third-order Hamiltonian 



1. Expansion of the condensate number operator 



III. CONSTRUCTION OF A THIRD-ORDER 
HAMILTONIAN 



A. Transformation of the full Hamiltonian 



Until now, every term with an explicit time-dependence has 
been shown with a t argument. From now on we neglect 
this, but it should be remembered that <p(r), A(r), N c , iV c , and 
<2(r, r'), all are in general explicitly time-dependent. 

One can readily transform the many-body Hamiltonian of 
Eq. (Q]), by everywhere expanding the field operators accord- 
ing to Eq. (O, and then collecting terms to produce products 



If the system is in a number eigenstate of total particle num- 
ber N, the number fluctuations of the condensate and non- 
condensate components must be equal and opposite. For- 
mally, 



(r) 



(30) 



= J drl^(r)^A(r)\- j drA + (r)^pA(r). 

To quadratic order in A(r), 

Nc=N c + j dr{A\r)A(r))- J drA f (r)A(r) (31) 

(the first corrections beyond this appear at quartic order and 
are not considered in this paper). To zeroth order N c = N c - 
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We now apply Eq. OTT ) to Eq. d29l ), keeping only terms of 
up to third order in the fluctuation terms. Pragmatically, this is 
equivalent to immediately abandoning the fourth-order term 
in Eq. d29j >, substituting N c for N c in the second- and third- 
order terms, and substituting Eq. (l3TT i into the zeroth- and first- 
order terms. This then produces: 



H 3 =N C drf(r) 



tfsp(r) + ||0(r)| 5 



<Kr) 



+ J dr{<f>*(r)[H sp (r) + C/|0(r)| 2 ] A(r) + H.c.} 
+ J drA\r) [// S p(r) + 2C/|^(r)| 2 ]A(r) 
+ | J dr[f(r) 2 A(r) 2 +U.c.]-j J dr\<p(r)f 
+ J dr' [<A + (r')A(r')> - A^rOACr')] 
x J" dr<p(r) [H sp (r) + U\<f>(r)\ 2 ] tfr) 
j dr [0*(r)A t (r)A(r) 2 + H.c] 
j dr[0*(r)|0(r)| 2 A(r) + H.c] 











ff drdr'{<p' 



(r)|0(r)| 2 



W c J J 

x [<A + (r')A(r')> - A t (r')A(r')] A(r) + H.c.}, 



(32) 



where the terms have been arranged in descending order of 
powers of V^Ve- 



An appropriate approximate Hamiltonian Hj consistent 
with this desired level of approximation in the equations of 
motion, should thus be such that the commutator [A(r), H3] 
produces terms contributing to the equation of motion for 
A(r) in the desired form. This means either scalar terms to 
zeroth order in the fluctuation operators, first-order operator- 
valued terms, or scalar second-order terms in the form of pair- 
averages. From Eq. ( l27l i and Eq. ( t3Tb it can be seen that, to 
quadratic order, 



[A(r), A'(r')] «fi(r,r')U+ dr 



,<At(r")A(r")> 



- dr 



„A T (r")A(r") 

Nr 



N c 

A t (r')A(r) 



(34) 



and that the first corrections appear at quartic order. To a 
Gaussian level of approximation, the operator products are 
consistently replaced by expectation values. Thus, 



[A(r), A 1 (r')] ~ Q(r, r ) . (35) 



N c 



When considering quadratic and cubic terms in the postulated 
third-order Hamiltonian H3, this commutator is simplified fur- 
ther to 



[A(r), A'(r')] = Q(r,r'), 



(36) 



2. Gaussian approximation of the fluctuation terms 

In the present context, a Gaussian approximation means 
that all expectation values of products of the fluctuation op- 
erators A(r), A' (r) are either zero (for odd products), or ex- 
pressable in terms of products of pair-averages OlOlll . Prag- 
matically, in the equation of motion derived for A(r), which 
we determine up to quadratic order in the fluctuation opera- 
tors, all quadratic products of A(r) and A ' (r) must be replaced 
by their expectation values. Doing this guarantees, for exam- 
ple, that a consistent Gaussian approximant to the equation of 
motion for the pair-average (A 1 (r)A(r')} is deduced directly 
from 



-<At(r)A(r')> 
dt 



d 



A(r') 



+ A"'(r) 



d - , 
dt Mr) 



(33) 



without there being any subsequent need for expansion of ex- 
pectation values of higher-order products of the fluctuation 
operators in terms of pair-averages [93]. 



as otherwise cubic and quartic terms appear in the final equa- 
tion of motion. 

Hence, we deduce that the cubic fluctuation operator prod- 
ucts appearing in H3 [Eq. ( l32b l must be expanded into sums of 
linear operator-valued terms multiplied by pair-averages. To 
this degree of approximation, this is accomplished by express- 
ing cubic products as the sum of all possible pair-averages, 
multiplied by the remaining fluctuation operator. This is 
equivalent to a Hartree-Fock factorization, as described, for 
example, in Ref. lf78tl . 

For example 



A" r (r)A(r')A(r") ^<A" r (r)A(r'))A(r") + <A T (r)A(r")>A(r') 



+ <A(r')A(r")>A 1 (r), 



(37) 



and we deduce that factorising the cubic terms appearing in 
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Eq. d32b results in: 



H 3 =N C drf(r) 



ff sp (r) + ||0(r)| 5 



<Kr) 



+ J rfr{0*(r)[// sp (r) + C/|0(r)| 2 ]A(r) + H.c.} 

+ J drA+Cr) [H sp (r) + 2U\cf>(r)\ 2 ] Mr) 

+ | J rfr[<f(r) 2 A(r) 2 +H.c.]--| J A#(r)| 4 

+ J dr' [<A + (r')A(r')> - A^OACr')] 

x J ^r0*(r) [// sp (r) + C/|0(r)| 2 ] 0(r) 



U 



j dr[q 



'(r) 



x [2<A t (r)A(r))A(r) + AVXACr) 2 )] + H.c.} 
--jj^f dr [<f (r)l<Kr)| 2 A(r) +H.c] 

+ -j=ff drdr'[f(r)\<f>(r)\ 2 

x [<A t (r')A(r))A(r') + A + (r')<A(r')A(r)>] + H.c.}. 



(38) 



It is with respect to this third-order Hamiltonian that our 
second-order equations of motion will be defined. 

The factorization procedure is analogous to that used in 
Hartree-Fock-Bogoliubov methods. As such, it is not gen- 
erally valid; careful consideration reveals this not to be a seri- 
ous problem in the present specific context, however. Hartree- 
Fock-Bogoliubov factorizations have also been applied in the 
full binary interaction Hamiltonian to both cubic and quartic 
products of the fluctuation operators. Work by Morgan ll57ll 
revealed that factorization of the cubic products omitted terms 
which were as large as terms of quartic origin which were re- 
tained. We, however, have already eliminated quartic fluctua- 
tion terms from consideration, and in the steady state all cubic 
terms will also be eliminated (see Section HVD 4b . If exten- 
sion of the theory to include higher-order terms is desired, this 
simplification will need to be revisited. 



IV. EQUATIONS OF MOTION 
A. General properties of the equations of motion 

1. Explicit time dependences 

It is convenient to have expressions describing the explicit 
time-dependence only of a c and 8*¥(r). 

Taking the partial time-derivative of Eq. ([Tol l, and re- 
calling that the bosonic field operator has no explicit time- 



dependence, we deduce that 

= J dr¥\r) 



d 

ifi—Mr) 
dt 



(39) 



Similarly, taking the partial time-derivative of Eq. (fTTT i pro- 
duces 



at 



J" 



m-Q(r,r') 
at 



<F(r'). 



(40) 



The condensate mode-function <f>(r) is defined to have unit 
norm, which directly implies 



dr 



dt 



4>*(r) 



(r) 



dt 



<P(r) 



(41) 



The resulting Eq. fiTb can then be substituted into Eq. ( l40b . 
producing 



i%— $$(r) = - a c f dr'Q(r,r' 
at J 



m—4>(r ) 

at 



4>{r) 



ih-J>*(r') 
at 



(42) 



di'ir'). 



In Eq. (|39l and Eq. ( l42l i. we have the final forms of the 
desired expressions. 



2. Condensate number 

The general equation of motion for the condensate number 
operator, N c - a\a c , is 



ih— = [N c ,H] + iti—, 
dt at 



(43) 



from which the dynamics of N c = (N c ) are deduced by taking 
the expectation value. 

We first consider the explicit time dependence in isola- 
tion. Substituting Eq. d39l and its Hermitian conjugate into 
dNddt = (da c /dt)a c + a c (da c /dt) produces 



0(r) 



dN c C ■> . n d 

to-gf = J dr[N c f(r)+ y/N c A\T)\ih- 

drift- f{r) [N c <p(r) + VA^A(r)] 

Substituting in Eq. ( |4TT i. we simplify Eq. d44b to 

M ... f 



(44) 



i%- 



dt 



N c I drA T (r) 
dr 



d 

ih—d){T) 
at 



ih-<l>*(r) 
at 



(45) 



A(r). 



Equation (05]) is entirely composed of linear fluctuation oper- 
ator terms. Hence, there is no explicit time dependence to the 
condensate number, i.e., 



dN c I dN c 

ih = ( ih 

dt \ dt 



0. 



(46) 



9 



Therefore, to all orders, the entire time-dependence of the 
condensate number follows from the (implicit) commutator 
term of Eq. ( |43T >: 



at 



(47) 



In principle this can be determined directly from the appropri- 
ate form of the Hamiltonian H. If one is in any case determin- 
ing the time-evolution of the individual fluctuation operators 
A(r), A 1 (r), it is generally more convenient to note from Eq. 
(T3Tb that N c - N - J<ir(A t (r)A(r)) to quadratic order, and 
therefore that 



th 



dN c = r 
dt " J 



dr ( A T (r) 



d ~ 
ih-Mr) 
dt 



d -, 
m— A'(r) 
dt 



A(r) , 



(48) 



to the (quadratic) order considered here. 



3. Fluctuation operator 



operator A(r), replace second-order terms with their expecta- 
tion values, and neglect higher-order terms altogether. Equa- 
tion ( IBTl l then simplifies to 



i%— Mr) = - 

dt 



dr'Q(r,r') 



«»-0(rO 
at 



dr' 



■ j dr'f(r') 



VAT 



dr' 



at 
d 

ifi—cf>(r') 



A(r') 



A(r) 



(52) 



<A 1 '(r')A(r)). 



Which of the terms of Eq. ( 1521 are subsequently retained 
depends on the order to which one wishes to carry out a given 
calculation. In order to determine the full dynamics to the 
desired order, we need to know the form of the appropriate 
approximate Hamiltonian. Sections [IV Bl IIV Cl and II V Dl de- 
duce such Hamiltonians to first, second, and third order, re- 
spectively, as well as the associated time-evolutions implied 
by them. 



We now consider the dynamics of the number-conserving 
fluctuation operator A(r) directly. In general, the Heisenberg 
time-evolution of the fluctuation operator is given by 



1 p\ 

m—k(r) = [A(r),H] + i%—k{r). 
dt at 



(49) 



We again initially consider the explicit time-dependence of 
Eq. d49l . which, from the definition of the fluctuation operator 
given by Eq. d25l ). yields 

8 8N 1 * - 

i%— A(r) = - ih—Z ald^ir) 

dt dt 2N C VA£ 



1 dal f 



(50) 



As there is no explicit time-dependence to N c [Eq. (|46]>1. the 
first line of Eq. (150b can be eliminated. After substituting in 
Eq. d39l and Eq. (l42l . what remains can be expanded in terms 
of fluctuation and condensate-number operators: 



^|-A(r) 
at 



-<P(r) J dr' 
dr'f(r') 



fi(r,r') 



d 

ih-J>(r') 
at 



d 

iti-f(r') 
at 

d 

ifi-<t>(r') 
at 

m^cp(r') 
at 



A(r') 
A(r) 



~ + N r ~ 

A^rO-^Afr). 

Nr 



(51) 



Working within the Gaussian approximation described in 
Section lTHB 21 we retain terms to first-order in the fluctuation 



B. Reduced first-order Hamiltonian 

1. Reduction to a first-order Hamiltonian 

In principle, one can consider a zeroth-order approximation 
to the Hamiltonian of Eq. (l38l . This is obtained by neglecting 
all fluctuation terms, and yields a classical energy functional 



H = N c J drf(r) 



tfsp(r) + ||0(r)| 2 



m- (53) 



The lowest order Hamiltonian of real interest to us is linear in 
the fluctuation operators A(r), A(r), which is when it first has 
a definite operator character. Dropping all terms of second and 
third order in the fluctuation operators from Eq. fl38l ) leaves the 
appropriate first-order form of the Hamiltonian: 



H x =N C drf(r) 



tf S p(r) + ||0(r)p 



<P(r) 



+ VaT j dr{<f>*(r)[H sp (r) + U\<p(r)\ 2 ]Mr) + n.c}. 

(54) 

2. Deduction of the Gross-Pitaevskii equation 

As we are using a first-order approximate Hamiltonian to 
deduce a zeroth-order approximate equation of motion, we 
combine Eq. (149b with the first line of Eq. (l52l [the other terms 
are neglected as being of linear or greater order in A(r)] , yield- 
ing 



ihj t kf) = [A(r),#i] " VA^c J dr'Q(r,r' 



(55) 
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Using the zeroth-order form of the commutator [Eq. (l36ill. 
inserting the first-order Hamiltonian [Eq. d54"l )l into Eq. d55l l 
produces 



ih— A(r) = Vw c I rfr'Q(r,r') 
of J 



H sp (r') + J/|^(r')| 2 - 



(56) 



0(r'). 



Taking the expectation value of Eq. 115 6b. and using the fact 
that (dA(r)/dt) = d{A(r))/dt = 0, we get the time-dependent 
Gross-Pitaevskii equation, in essentially the same manner as 
Castin and Dum 14511 . with N c taking the place of N (U = 

U N C ): 



ik^m = K>M + U\4>(r)\ 2 ~ *o] flr), (57) 



where 



i/ sp (r) + f/|0(r)| 2 -^- 



0(r). (58) 



By norm conservation [Eq. d4"Tl)l. the scalar value Aq = AX, 
and is therefore always real. Substituting Eq. d57| i into Eq. 
d56b then directly implies that ifid A(r) I dt = 0, and hence 
[through Eq. (|48]l] that dN c /dt = 0, i.e., to this order, there 
is no time-dependence to the non-condensate component, and 
hence no change in the number of non-condensate atoms. 

This is to be expected, as we are considering the system 
dynamics to zeroth order in the fluctuation operators. Thus, 
to the current order, we are ignoring the fluctuation operators 
altogether in the equations of motion. 



3. Time-independent case 

Assuming <p(r) to be a steady state with respect to Eq. (|57T i 
(generally, although not necessarily the lowest energy steady 
state), one derives the time-independent Gross-Pitaevskii 
equation 



Ao<p(r)= [tf sp (r) + W(r)| 2 ]0(r), 



(59) 



where Aq takes the form of a nonlinear eigenvalue, which at 
this level of approximation can be identified with the chemical 
potential. A consequence of this is that the linear terms in the 
first-order Hamiltonian [Eq. d56l ll can be eliminated, reducing 
H\ to the zeroth-order form given in Eq. 



C. Reduced second-order Hamiltonian 



I. Reduction to a second-order Hamiltonian 



Dropping all terms cubic in the fluctuation operators, A(r) 
and A T (r), from Eq. d38l > yields 



#2 



=N C J drf(r) 



tfsp(r) + y |#r) 



4>(r) 



dr\<f>(r)f 



+ ff ^r'<A t (r')A(r'))0 t (r)[// sp (r) + f/|^(r)| 2 ]^(r) 

+ j dr{<f>*(r)[H sp (r) + U\<p(r)\ 2 ]Mr) + H-C.} 

+ J drkUr) [// sp (r) + 2f/|^(r)| 2 ]A(r) 

+ jj dr[0*(r) 2 A(r) 2 +H.c] 

- Jj drdr'A\r')A(r')f(r) [# sp (r) + U\cf>(r)\ 2 ] «r), 



(60) 



where the terms have been arranged such that all scalar terms 
come first (including fluctuation operator pair-averages), fol- 
lowed by terms linear in the fluctuation operators, and subse- 
quently by quadratic (non-expectation value) fluctuation op- 
erator terms. 



2. Deduction of the modified Bogoliubov-de Gennes equations 



To determine the equation of motion for the number- 
conserving fluctuation operator A(r) to linear order, we must 
include the zeroth- and linear-order terms from Eq. d52l >. in- 
serting these and the quadratic Hamiltonian [Eq. (l60l il into Eq. 
(US: 



^A(r)=[A(r),£ 2 ]- ^JW C J dr'Q(r,r' 



m J dr' 
J dr'f(r' 



d 

m-J>*(r') 
ot 

d , 
m—(p(r ) 
ot 



A(r') 
A(r). 



(61) 



We continue to use the zeroth-order form of the commutator 
[Eq. 06H . as to this order we may still neglect the quadratic 
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correction. Applying this to Eq. doTT ) yields 

ifi— A(r)=V# c f rfr'Q(r,r') 
a? J 



ff sp (r') + f7|0(r')| 2 - ift- 



(O 



J fl fr'e(r,r')[i/ sp (r') + 2j/|^(r')| 2 ]A(r') 
J rfr'e(r,r')A t (r')0(r') 2 

J dr'f(r' 



+ U 
~<P(r) 
-Mr) 



iH—cf>*(r') 



A(r') 



H sp (r') + U\4>(r')\ 2 - ihj t 



(r'). 



(62) 



Taking the expectation value produces the same Gross- 
Pitaevskii equation [Eq. d57b l deduced in Section lTV B 21 This 
is due to the fact that no linear terms not already present in 
Eq. ( |54b appear in Eq. d60l >. 

Equation ( T57b can be substituted back into Eq. i62i , simpli- 
fying it to 



^-A(r) = [H sp (r) + U\<p(r)\ 2 ~ 4>] Mr) 



dt 



+ U 

+ u 



J dr'Q(r,r')\c/>(r')\ 2 A(r') 
J dr'Q(r,r')<f> 2 (r')A\r'). 



(63) 



Equation ( 1631 . together with its Hermitian conjugate, form the 
Bogoliubov-de Gennes equations ll4ll l42ll . modified slightly 
by the presence of the orthogonal projectors Q(r, r'). This is 
equivalent to the result presented by Gardiner [44] and Castin 
and Dum |45 |. apart from the use of N c rather than N. 

The presence of the projectors is due to the fact that the 
definition of the condensate and non-condensate components 
[Eq. (O] explicitly guarantees their orthogonality 114311 . This 
is not true with a conventional symmetry-breaking approach. 
Note, however, that if one considers a spatially homogeneous 
condensate density, then 

i*j t Hr) = K(r) + 2U\<p(r)\ 2 ~ <*o] Mr) + U<p 2 (r)M(r). 

(64) 

which coincides with the conventional form of the 
Bogoliubov-de Gennes equations l4ll l42ll . 

3. Number evolution 



Substituting Eq. ( I63t . together with its Hermitian conju- 
gate, into Eq. (l48l yields that the condensate number evolves 
as 

i% lt = ° I dT ^* (r)2<A(r)2> - < A "' ( r ) 2 >^( r ) 2 ] • (65) 



This equation is composed of terms quadratic in the fluc- 
tuation operators, even though we have everywhere else ne- 
glected equivalent quadratic terms. One could argue that 
these contributions should be consistently neglected as being 
"small" compared to the current (linear) order of interest. In- 
consistencies, such as the fact that when considering a non- 
steady-state evolution the modified Bogoliubov-de Gennes 
equations [Eq. (l63l l can imply unconstrained growth of the 
non-condensate fraction without there being any correspond- 
ing effect on the condensate evolution [Eq. (fSTj ll ll26ll27[|28ll . 
would then be dismissed as being due to a misguided attempt 
to use a first-order theory to determine the evolution of a 
second-order quantity. It is useful to examine what is going 
on in a little more detail, however. 

Caveats about purely quantum fluctuations aside, finite 
fluctuations generally imply finite pair averages, as discussed 
in Section III C 1 1 and Section III C 41 Fluctuations assert their 
existence by having an observable effect, through pair aver- 
ages (as well as possibly higher-order connected correlation 
functions) [100]; for example, when determining the quan- 
tum depletion 0,181]. For a fixed, finite total number of parti- 
cles there is, therefore, a consistency problem intrinsic to first- 
order theories in the fluctuation operators. In Section ITV Dl it 
will be shown that, while extending the theory to second-order 
modifies the equation of motion for <p{r), it does not change 
the form of the equations of motion for A(r) [Eq. (l63ll. and 
hence A^ c [Eq. d65l>l. As will be described in Section irVD5l in 
an infinite particle limit (N c — > oo while Af, remains finite) the 
first-order equations of motion of Section [TV Cl are recovered 
from the second-order equations of motion. Strictly speaking, 
this limit is required for the first-order equations of motion to 
be meaningful. 

When allowing N c (and hence N) to tend to infinity while 
requiring N, to remain finite, it is preferable to consider the 
time evolutions of Af, and N c /N, which are both finite. Instead 
of Eq. d65l ), we then have 



iti— (—)=U f dr 
dt\N I J 



t(r)2 ^_<^f) 0(r)2 



N 



N 



(66) 



which tends to zero on the right-hand-side as — > oo, and 



dN, 
dt 



U 



j dr [<A + (r) 2 >0(r) 2 - 0*(r) 2 <A(r) 2 > 



(67) 



which does not. Within this limit it is perfectly legitimate to 
use Eq. (|67| | to determine the time-evolution of the number 
of non-condensate particles; change in N, then corresponds to 
an infinitesimal fractional change in A^ c . Conversely, Af c /Af = 
1 - N,/N — > 1, and so the time-derivative in Eq. d66l ) must 
logically be zero. 

Such a limit is not appropriate for the desired self-consistent 
propagation of condensate and non-condensate dynamics with 
a fixed, finite total particle number N, for which the full 
second-order treatment described in Section |IV D| is neces- 
sary. 
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4. Time-independent case 



2. Deduction of the generalized Gross-Pitaevskii equation 



As in Section HV B 31 we substitute the time-independent 
Gross-Pitaevskii equation [Eq. d59l )l into the second-order 
Hamiltonian [Eq. d60ill, and eliminate the same linear terms. 
This yields a form of the Hamiltonian, 



H 2 



*(r) 



tfsp(r) + ||0(r)p 



0(r) 



(68) 



=N C j dr<f>*( 
+ A Q j dr(A l (r)A(r)} - | J A#(r)| 4 
+ J rfrA^r) [ff S p(r) + 2f7|0(r)| 2 - Ao] A(r) 

+ dr[^(r) 2 A(r) 2 +H.c.], 

equivalent to that deduced in a number-conserving fashion by 
Gardiner JH. 



D. Properties of the third-order Hamiltonian 

1. Gaussian form of the third-order Hamiltonian 

The appropriate Gaussian third-order form of the Hamil- 
tonian is exactly as given in Eq. ( 1381 ). As in Eq. d60l l. it is 
convenient to rearrange the equation such that all scalar terms 
come first (including fluctuation operator pair-averages), fol- 
lowed by terms linear in the fluctuation operators (including 
those multiplied by fluctuation operator pair-averages), and 
subsequently by quadratic fluctuation operator terms: 



H 3 =N C drf(r) 



tfs P (r)+^Wr)f 



0(r). 



j A#(r)| 4 

+ ff rfrrfr'(A t (r')A(r')>0*(r)[H sp (r) + C/|0(r)| 2 ]0(r) 
+ j dr [tf(r) [tf sp (r) + U\<f>(r)\ 2 ] A(r) + H.c.} 

u r 

+ —= dr f20*(r)<A 1 "(r)A(r)>A(r) + H.c.l 
+ ^j=jdr [<A 1 '(r) 2 )0(r)A(r) + H.c] 

--j=ff drdr ' {l^( r )l 2 [0*(r)<A t (r')A(r)> 
+<A t (r')A t (r))0(r)]A(r') + H.c.} 
-^j=jdr [<f (r)|0(r)| 2 A(r) + H.c] 

+ J drkUr) [i/ sp (r) + 2f7|^(r)| 2 ]A(r) 

+ ~ J rfr[0*(r) 2 A(r) 2 +H.c] 

- ff rfr c /r'A t (r')A(r')0*(r) [// sp (r) + U\<p(r)\ 2 ] 

(69) 



We now determine the equation of motion for the number- 
conserving fluctuation operator A(r), to quadratic order. We 
substitute Eq. d52l ). in its entirety, and the Gaussian form of 
the cubic Hamiltonian [Eq. d69l )l into Eq. j49l . The equation 
of motion can then be written as: 



m^-Mr) =[A(r), 
dt 



H 3 ]~ V^Vc J dr' 

-0(r) J dr' 



d 
ot 



j dr'f{r') 



d 



<A t (r')A(r)) 



ihj-<t>*{r') 
ot 



A(r') 



A(r). 



(70) 



To produce a consistent second-order equation of motion, we 
must now include the quadratic correction to the fluctuation 
operator commutator, using the full form given by Eq. $35[ . 
This will also produce cubic and quartic corrections, which 
should be consistently neglected. Effectively this means that 
we use the full form of the commutator when determining the 
time-dependence due to terms of Eq. d69l that are linear in the 
fluctuation operators. Otherwise, the zeroth-order form [Eq. 
(|36| >1 will suffice. 

Doing this produces, subsequent to some rearrangement, 



d ~ 
zfc— A(r) 
dt 



■^NcJ dr'{ 



g(r,r') // S p(r') 



+ U 



Va c 



<A t (r')A(r)> 



H sp (r') + 2f7|0(r')| 2 - ih 



<f>(r') 



+ f/0*(r')<A(r')A(r))|^(r')| 2 1 

+ j dr'Q(r,r'){[H sp (r') 
+ 2f7|0(r')| 2 ]A(r') + Uk\r')<f>(r') 2 } 
A(r') 



(71) 



-0(r) j dr' 
-Mr) j dr'fir' 



iti-J>*(r') 
ot 



H sp (r') + U\<f>(r')\ 2 - ih- 



(r')- 



As in Section ITVC 21 taking the expectation value of this ex- 
pression eliminates all the linear fluctuation terms, leaving us 
with an equation of motion for the condensate mode <p(r). Un- 
like the simple Gross-Pitaevskii equation [Eq. (|57|)1. this equa- 
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tion of motion couples to normal and anomalous pair-averages 



Note that A 2 , unlike Aq [Eq. d58l>l. may be complex. The first 
integral, in a similar fashion to Aq, can be seen to be always 
real. This is not necessarily so for the second integral, as can 
be seen from 

A 2 - XI = ±-U fdr [«f (r) 2 <A(r) 2 > - (A^r) 2 )^) 2 ] . (74) 



the final form of the generalized Gross Pitaevskii equation 
(formally, in conjunction with its complex conjugate, to which 
it is coupled). This is essentially as was used by Morgan J39tl 
to explain finite temperature effects on the excitation spec- 
trum measu red i n the JILA 87 Rb Bose-Einstein condensate 
experiment Il02ll . The anomalous average (A(r) 2 ) must be 
appropr iately renormalized to avoid ultraviolet divergences 
J53,|6ll|2i|9l|9l|97l|9itl,asis briefly sketched in Appendix 
El 

In Ref. i45ll . Castin and Dum choose to describe the con- 
densate mode in terms proportional to inverse powers of the 
total particle number, i.e., as 

</>(r) = o (r) + —0i (r) + ±fe(r) + ■ ■ • . (76) 

n 

The zeroth-order term (f>o(r) is then propagated by the Gross- 
Pitaevskii equation [Eq. (|57T >. but Castin and Dum formally 
have U = UqN c replaced by UqN]. The first-order term (p\(r) 
is equal to zero [as described in Section IIVC21 there are no 
first-order corrections to the time evolution of </>(r)]. Equa- 
tions (95) and (96) of Ref. ll45tl describe the time-evolution 
of the component of 02(f) orthogonal to </>o(r), which are cou- 
pled to the motion of (f>o(r) propagated by the Gross-Pitaevskii 



of the number-conserving fluctuation operators: 



(72) 



(73) 



We can eliminate the time-derivative on the right-hand side 
of Eq. (1721 by iterative resubstitution, keeping only terms of 
up to the appropriate order. This is equivalent to substituting 
in the lower-order equation of motion for <p(r), i.e., the Gross- 
Pitaevskii equation [Eq. (|57|>1. Doing this produces 



(75) 



equation. Combining these equations appears to be equivalent 
to substituting <p(r) = (f>o(r) + </>2(r)/N into Eq. d75l l, dropping 
terms in (f> 2 (r) considered to be of too high an order, and re- 
placing 1 /N c with 1 IN. We have found it simpler to consider 
the evolution of a unified <p{r) directly, and the deduction of 
Eq. ( |75l l is a significant step toward the paper's main result. 

Noting that J dr'Q(r, r')<A 1 "(r")A(r')) = <A f (r")A(r)> and 
that similarly J"rfr'g(r,r')<A(r")A(r')> = <A(r")A(r)>, we 
see that substituting Eq. ( l72l into Eq. ( |7T1 i. the equation of 
motion for A(r), causes all terms not linear in the fluctua- 
tion operators to vanish. This is basically equivalent to the 
removal of the zeroth-order terms when deducing the modi- 
fied Bogoliubov-de Gennes equations in Section lrVC2l 

One can again substitute Eq. d75l l for ihd<p(r)ldt where it 
appears in what remains of Eq. dTTl l, neglecting all higher or- 
der terms; note, however, that this is equivalent to substituting 
in the simple Gross Pitaevskii equation [Eq. d57lil. This leaves 
us with the same modified Bogoliubov-de Gennes equations 
[Eq. d63l )l as determined previously, in Section [rVC2| 

The generalized Gross-Pitaevskii equation [Eq. (|75|>1. to- 
gether with the modified Bogoliubov-de Gennes equations 
[Eq. j63l l thus describe the second-order coupled condensate 
and non-condensate dynamics, respectively. This states the 



ih-(KT)=\H tp (r) + u 



1-— wmf+2 



<A'(r)A(r)> 



A 2 \ 4>(r) + Uf(r) 



(Mr) 2 } 



<At(r')A(r)) 

N c 



H s Jr') + 2f/|0(r')| 2 - ih 



, >s fiM*< ^ |2 (A(rQA(r)) 
(r') + U<f> (r )|0(r')r t; 

l*c 



where the scalar value, A 2 , is given by 

A 2 = J drf(r)(H sv (r) + 



1 \ , (A'(r)A(r)) 

1 \ m ? + 2L-yi» 



-ih-}4>(r) + U 



f, ... , 2 <A(r) 2 ) 
J dT(t> (r) " — ' 



r 



1 * r 



ifi-m = I ff sp (r) + U ||1 - — I |0Cr)r + 2 



<A t (r)A(r)> 



Nr. 



A 2 0(r) + U<f>*(r) 



(Mr) 2 } 

Nc 



-0 j dr'|0(rO 



<At(r')A(r)> <A(r')A(r)> 
■<p(r ) + 4> (r ) 



Nc 



N r 
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main result of the paper. 

It should be emphasized that the evolution predicted by the 
modified Bogoliubov-de Gennes equations may be very differ- 
ent if it is coupled to the generalized Gross-Pitaevskii equa- 
tion [Eq. (l75l) l rather than the simple Gross-Pitaevskii equa- 
tion [Eq. (|57|>1. That this may constitute a more consistent 
treatment is shown by the fact that, just as there is an action 
of the condensate normal and anomalous density terms on the 
time-evolution of the number conserving fluctuation operators 



[Eq. (l63lll. there is a corresponding back action of the normal 
and anomalous pair-averages on the time-evolution of the con- 
densate mode [Eq. (l75lll. 

A similar generalized Gross-Pitaevskii equation can be de- 
rived within a symmetry -breaking context 117 811 . but without 
the integral term on the second line of Eq. d75l l. Before 
discussing the role of this term, we note that the projectors 
Q(r, r') in the modified Bogoliubov-de Gennes equations [Eq. 
(l63l ll can be expanded to give 



foj t kr) = K(r) + 2U\cf>(r)\ 2 ~ 4>] Mr) - tfrftfOr) ~uj dr'\4>(r')\ 2 [<f (r')0(r)A(r') + A'(r')0(r')0(r)] . (77) 



r 



Those parts of the integral terms of Eq. (l75l l and Eq. tlli en- 
closed within square brackets are of almost identical form, but 
with the roles of the condensate mode functions and the fluc- 
tuation operators exchanged. A comparably elegant simplifi- 
cation of notation afforded by use of the projectors in Eq. (l63l l 
is not obvious for Eq. d75l l. The function of the integral terms 
in Eq. ( fTTT i and Eq. d75l l is equivalent, however — to ensure 
that the orthogonality of the condensate and non-condensate 
components is maintained j45|] . Hence, their explicitly nonlo- 
cal form, and the consequence that both integral terms vanish 
in the limit of a spatially homogeneous condensate density. 

The appearance of such a term at this order is necessar- 
ily in conjunction with the coupling of the generalized Gross- 
Pitaevskii equation [Eq. d75l )l to the fluctuation operator nor- 
mal and anomalous densities. This is unlike the case in 
Section IIV CI where the result of the simple time-dependent 
Gross-Pitaevskii equation [Eq. (l57l il feeds into the modified 
Bogoliubov-de Gennes equations [Eq. d63lll. but the Gross- 
Pitaevskii equation itself evolves in complete isolation. 



3. Number evolution 

As the time-evolution of the number-conserving fluctua- 
tion operators, A(r) and A' (r), is still given by the modified 
Bogoliubov-de Gennes equations [Eq. d63lll. the condensate- 
number evolution must still be given by Eq. (f65t . Note, how- 
ever, from Eq. (l74l i. that the number dynamics can also be cast 
as 



dN c 
dt 



Ai - X% 

ih 



l -N e . 



(78) 



This has the form of a simple linear differential equation. The 
(time-dependent) rate of growth or decay of the number of 
condensate particles is equal to the difference between the cre- 
ation of pairs of condensate particles in conjunction with the 
annihilation of pairs of non-condensate particles, and the re- 
verse process. 

The significance of this result is that the condensate-number 
evolution equation directly implied by the third-order Hamil- 
tonian contains no terms of higher than second-order in the 



number-conserving fluctuation operators, which is consistent 
with the order of those fluctuation-operator terms appearing in 
the generalized Gross-Pitaevskii equation. This is the lowest 
non-trivial order at which such a consistent description is pos- 
sible for a finite number of particles ifiol . One might have 
expected higher-order fluctuation operator terms to be neces- 
sary in the non-condensate evolution for a treatment consis- 
tent with the generalized Gross-Pitaevskii equation ll75l . This 
is not so; consistent number dynamics in fact require that there 
be no extension to the modified Bogoliubov-De Gennes equa- 
tions [Eq. d63lll. 



4. Time-independent case 

If we assume a steady state for 0(r), then the [equivalent 
to Eq. d59l ll time-independent generalized Gross-Pitaevskii 
equation is given by 



A 2 <p(r) =\H sp (r) + U 

,<A+(r)A(r)> 



+ 2- 



I 



- dr' 



<A'(r')A(r)> - 

■U\cf>{r')\-4>(r) 



(79) 



+ f(r')U\<p(r')\ 



, ,<A(r')A(r)> 



Substituting this back into Eq. j69l , all linear and cubic-order 
terms disappear. This is analogous to the way all linear-order 
terms disappeared in the derivation of the second-order time- 
independent Hamiltonian, and the elimination of these terms 
leaves us with the same form of time-independent Hamilto- 
nian [Eq. d68ll. 



5. Infinite-particle limit 

Examination of Eq. (p75T > and Eq. d77| i reveals that allowing 
the number of condensate particles to arbitrarily increase, i.e., 
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N c — > oo, causes all higher-order terms present in the gen- 
eralized Gross-Pitaevskii equation to vanish, leaving the sim- 
ple Gross-Pi taevskii equation [Eq. (l57lil. whereas the modified 
Bogoliubov-de Gennes equations are unchanged. 

We thus reduce exactly to the first-order formulae gained 
using an approximate second-order Hamiltonian [Eq. <f60b 1 . 
When one considers that a treatment using the modified 
Bogoliubov-de Gennes equations [Eq. d63l )l coupled to the 
simple Gross-Pitaevskii equation [Eq. (l57b l allows for unlim- 
ited growth of the non-condensate fraction without there being 
any effect on the condensate dynamics, it is clear that only in 
the limit of an infinite number of condensate particles can the 
dynamics predicted by these equations be strictly correct. 



V. EQUILIBRIUM PROPERTIES 
A. Overview 



Section IV Bl recaps the situations described by Gardiner 
l44ll and Castin and Du m l45ll, which, in addition to work by 
Girardeau and Arnowitt 114a, [4711 . sought to provide a number- 
conserving equivalent to the symmetry breaking Bogoliubov 
formalism l4ll l42tl . That is, considering the Hamiltonian to 
second order in the fluctuation terms, or equivalently, equa- 
tions of motion of up to first order in the fluctuation terms. In 
the present context, this is equivalent to assuming the correct- 
ness of Eq. ( |57| | and Eq. ( 1631 . Having set context and notation, 
Section lVCl considers some of the difficulties in going beyond 
this level of approximation. 



As Jc/r'2(r,r')A(r') = A(r), we choose to describe the 
fluctuation operator time-evolution by 



it<-\ tr'. 1 = J dr'£(v. i-'ij tV.V, I- (84i 



d_t A(r) 
'dt\ A + (r) 



.>J Mr') 
At(r') 



where 



X (r,r') = | f£ r \ M X y \\, (85) 
1 -M (r,r ) —L (r, r ) / 



and the elements of £(r, r') are defined by 

L(r, r') =6(r - r') [# sp (r') + U\cf>(r')\ 2 - A ] 

+ J dr"Q(r,r")U\<f>(r")\ 2 Q(r",r'), 

M(r,r') = j dr"Q(r,r")U<p(r") 2 Q*(r",r'). 



(86) 
(87) 



Note that L(r',r) = L*(r, r'), i.e., L is Hermitian, and thus 
X(r, rQ has some symmetry properties which JT(r, r') does 
not 14511 . 

Inserting the projectors Q(r,r') into Eq. (f80b in this way 
has the useful property that the evolutions predicted by the 
modified Bogoliubov-de Gennes equations [Eq. d63l >l and the 
simple Gross-Pitaevskii equation [Eq. d57l il are unified by the 
application of the operator X(r, r') onto an appropriate two- 
component state. Thus, replacing (A(r'), A' (r')) in Eq. (l84l i 
with (<p(r), 0) or (0, <p*(r)) reduces it to the simple Gross- 
Pitaevskii equation, or its complex conjugate, respectively 



B. Quasiparticle formulation 

1. Two-component representation 

As the time-evolution of the number-conserving fluctuation 
operator A(r) [Eq. d63l )l causes it to couple to its Hermitian 
conjugate, it can be convenient to write the coupled time- 
evolution equations in a unified two-component form. Thus 



where 



J(r,r') 



J(r, r') K(r, r') 



-K*(r,r') -7*(r,r') 

and the elements of J'ir, r') are defined by 

7(r, r') =«5(r - r') [# sp (r') + U\cf>(r')\ 2 - A ] 

+ Q(r,r')U\cf>(r')\ 2 , 
K(r,r') =e(r,r')f/0(r') 2 , 

and their complex conjugates. 



(81) 



(82) 
(83) 



2. Quasiparticles 
The spectral decomposition of £.(r, r'), 

X(r ) r0=|je i (^)(4(r'),-vl(r')) 



(88) 



the derivation of which is outlined in Appendix iBl provides a 
useful basis in which to expand the number conserving fluc- 
tuation operators: 

In turn, using the orthonormality relations 

J dr[u* k ,(r)u k (r) - v^(r)v t (r)] =6 W , (90) 
J dr[u k ,(r)v k (r)-v k ,(r)u k (r)]=0, (91) 
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we determine that the operator coefficients are given by 

b k = J fl fr«*(r)A(r) - v^(r)A t (r) (92) 



b' k = | dru k (r)A 1 (r) - v*(r)A(r) ( 93 ) 

and that their commutation relations are 

[h,bl,]= jj drdr'[ul(r)u k ,(r')-vl(r')v k ,(r)] 



x[A(r),AV)L 
[bM = ff drdr'[ M *(r)v*(r')-v*(r')4(r)] 



(94) 



(95) 



X[A(r),AV)L 



If we can assume the commutator for the number- 
conserving fluctuation operators to be reduced to the projec- 
tor Q(r,r') [Eq. d36bl. the operator coefficients b k , b\ form a 
bosonic algebra: 



[bk, b\,] =6w 

[h,h>] =0. 



(96) 
(97) 



The operators b, and b k are then quasiparticle creation and 
annihilation operators 11441 14511 . 



3. Reformulation of the Hamiltonian in terms of quasiparticles 

Substituting Eq. d89l into Eq. d68l [and making use of Eq. 
( IB21 . Eq. <H3]), Eq. dH, Eq. ((M), and Eq. (H)] yields 

oo ~ 

fi 2 =<H+ J] (^^-Yblh, I dru k (ry Uk ,(r) 



where 
<H=N, 



oo _ 



ff sp (r)+|(l-i)| ( /»(i-'- 



(98) 



oo 

oo p 
oo ~ 

^0 X (b\bl) dru* k (r)vl, 

k,k'=l 



(r) 



<irv^(r)v^,(r) 
dru k (r)u kl {r) 
drv k (r)u k ,(r) 
(r). 



(99) 



Making use of Eq. (1961 1, i.e., assuming the quasiparticle oper- 
ators to have bosonic commutation relations, reduces Eq 



to diagonal form [41|,|44 

OO n OO 

di-'H — ^ e k I </r|v,(r)| 2 + £ (100) 
and assuming a thermal equilibrium state, 'H reduces to 



ff sp (r) + -|l~ — ||0(r)| 



0(r) 



=N C J~ dr<p*(r) 
+ /1 ^]<^) I <ir [ M *(r) Mt (r) + v^(r)vt(r)] ( 10 1) 
+ ^o£ f drvK 

fc_ 1 *J 



t(r)v k {r). 



This all being so, the quasiparticle populations for a system 
in thermal equilibrium are given by (b k b k ) = [exp({e^ - \p, — 
M)WlkbT) - l] -1 l4l l39ll . where fi is the chemical potential, T 
the temperature, and kg Boltzmann's constant. Having pop- 
ulated the system appropriately, one can determine the time- 
evolution of the fluctuation operators from a system initially 
at equilibrium purely through the mode functions, such that 

and the b k , b\ are constant. 



C. Further considerations 



As has been shown in Section IIV C 41 and Section IIV D 41 
ft?, and H2 have the same form if the system is in equilib- 
rium [Eq. (I68H . meaning that Eq. d98l is an equally valid re- 
formulation of H3 in an equilibrium context. A concern is 
that use of the more complete formulation of the commutator 
[Eq. d35l l reveals that the quasiparticle commutation relations 
are not exactly bosonic [Eq. d94l ) and Eq. ([95])]. If we recall 
that, in conjunction with second-order terms, we have always 
used the simpler form of the commutator, in the context of the 
present paper this does not seem to be a critical consideration. 
Extending this approach to a consistent higher-order formal- 
ism, as is in principle desirable, will present some difficulties, 
however. 

We could take slightly different operators, defined from 
A c (r), Aj(r) [Eq. C3], 



A c (r) 
AJ"(r) 



5> 

k=\ 



u k (r) 
v*(r) 



SI 

k=\ 



(103) 



As a consequence of the commutation relation described in 
Eq. j23l , the commutation relations of b k and b k are exactly 
bosonic, and therefore could potentially better describe the 
system in terms of a Bose-Einstein distribution. This would 
be more in keeping with the spirit of the detailed treatment, 
making use of second-order perturbation theory, given in Ref. 

m. 
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As described in Section III CI A f (r) is not a simple fluctu- 
ation operator. In other words, its expectation value is not 
defined to be exactly equal to zero. The formal development 
of the equations of motion in Section [IV] explicitly relies on 
the number-conserving fluctuation operator having expecta- 
tion value equal to zero. Hence, using A c (r), which has ex- 
pectation value only approximately equal to zero, would mean 
that corrections, particularly at higher order, would have to be 
carefully calculated and accounted for. Furthermore, the idea 
that one is looking at fluctuations about a well defined mean, 
and that the magnitude of these fluctuations is observed, in 
the first instance, through pair expectation values, possibly to 
be followed by a well-defined hierarchy of connected correla- 
tion functions or cumulants, is muddied. In other words, an 
imprecision in the procedure used to determine the relevant 
equations of motion is introduced at a basic level. It therefore 
does not seem that a demand for such precision is compatible 
with defining perfectly bosonic quasiparticle operators. 

It should be emphasized that this is not a comment over 
whether it is more correct to use either A c (r) or A(r). Hav- 
ing an expectation value exactly equal to zero, and inducing 
exactly bosonic quasiparticle commutation relations, are both 
highly desirable properties for a hypothetical fluctuation oper- 
ator to have. Within the formal framework used in this paper, 
for a finite total number of particles, it seems that one has 
choose which of these properties is the more important for the 
purpose at hand. The emphasis in this paper is on the formal 
development of the relevant equations of motion. For this, at 
least, describing things in terms of A(r) and A'(r) seems to 
be more convenient. For other purposes, it may be more ap- 
propriate to consider a formulation closer to A c (r) and AJ(r) 

Good results have been achieved, using the equations of 
motion developed in this paper, in describing excitations in fi- 
nite temperature Bose-Einstein condensate by Morgan l39ll . 
We also note that such issues as potentially imperfecfly- 
bosonic quasiparticle operators are largely avoided if the ini- 
tial system has a negligible non-condensate fraction, even if 
subsequent dynamics (fo r ex amp le investigations of chaotic 
dynamics JM |H HI M, MM^lM, El El [M El) can 
cause significant depletion 112a 1271 12811 . hence requiring the 
kind of self-consistent treatment presented here. 



VI. CONCLUSIONS 

In conclusion, we have shown that a coupled system of 
equations, the generalized Gross-Pitaevskii equation and the 
modified Bogoliubov-de Gennes equation are the necessary 
minimally complete description to imply internally consis- 
tent number dynamics for a finite total number of particles. 
In other words, dynamics such that only particles lost from 
the condensate fraction are assumed by the non-condensate 
fraction, and vice-versa. Elaboration of the (linear) modified 
Bogoliubov-de Gennes equations is neither desirable nor nec- 
essary, as this would automatically lead to inconsistent num- 
ber dynamics. That an approach to second order in the fluctu- 
ation operators is necessary is directly implied by elementary 



statistical considerations; effectively that a finite fluctuation 
directly implies a finite variance, or its equivalent. Hence, in 
an infinite particle limit the first-order approach, consisting of 
the simple Gross-Pitaevskii equation coupled to the modified 
Bogoliubov-de Gennes equations, is recovered. It is only in 
this limit that the dynamics predicted by this system of equa- 
tions are technically consistent. A similar form of the ap- 
proach presented here has been employed 11381 l39l 14011 as a 
key component of an analysis of the observed excitations in 
finite temperature Bose-Einstein condensates, to good agree- 
ment with experiment j25ll . The formalism presented here 
is also suitable for the study of dynamically unstable Bose- 
Einstein condensate dynamics, where, even if the sample is 
initially at zero temperature, it is possible for a sizable non- 
condensate fraction to build up over time. 
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APPENDIX A: RENORMALIZATION OF THE 
ANOMALOUS AVERAGE 

The generalized Gross-Pitaevskii equation [Eq. ( |75T ll con- 
tains the anomalous average (A(r)A(r)), which is ultra-violet 
divergent. We give a brief summary of the reason and cure for 
this problem 03, HI El |H M M HI • 

The divergence arises from of the use of the contact poten- 
tial approximation. A genuinely ab initio theory would start 
by describing particle interactions using the true two-body po- 
tential. The contact "potential" is rather the zero-momentum 
limit of the two-body T-matrix describing the scattering of two 
particles in vacuum. It is introduced at the outset [Eq. (JT)] 
for a number of reasons: partly for convenience, partly be- 
cause this is the experimentally relevant quantity, and partly 
because it makes sense to include as much two-body physics 
as possible before embarking on a difficult many-body calcu- 
lation. We certainly cannot treat the two-body interaction with 
perturbation theory. This is apparent from the fact that the in- 
teractions can be described by a contact potential dependent 
only on a scattering length, whereas a perturbative treament 
would depend on the details of the potential [104]. 

However, this does mean that we have implicitly included 
at the outset various physical effects which must also appear in 
the many-body treatment. To avoid double-counting we need 
to subtract off the perturbative approximation to the two-body 
effects whenever we encounter them. 

The leading order interaction term is the nonlinear term in- 
volving the condensate. The interaction strength Uq in this 
expression must now be replaced by the second order approx- 
imation, i.e., the U in f/|0(r)| 2 0(r) must be replaced in Eq. 
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(75) by U + AU/N C , where 



AU 



u r « 



<m) 2 



(Al) 



and AU/N^ is the second order correction to the interaction 
strength as calculated from the Lippmann-Schwinger equa- 
tion. This correction can be grouped with the term in the 
generalized Gross-Pitaevskii equation [Eq. (175) 1 involving 
the anomalous average. This leads to a finite renormalized 
anomalous average m R {r), defined by 



m\r) = <A(r)A(r)> + -^-<p(r) z - 



AU 



(A2) 



It should therefore be implicitly assumed that the anoma- 
lous average (A(r)A(r)) appearing in Eq. (75) is replaced 
by m R (r) to produce a consistent, renormalized generalized 
Gross-Pitaevskii equation. 



APPENDIX B: SPECTRAL DECOMPOSITION OF £(r, r ) 

The treatment described in this appendix echoes that of 
Castin and Dum M45I1 . and is included for the sake of com- 
pleteness. 

We assume that (u k (r'), v k (r'J) is a right eigenstate of the 
operator £.(r, r') defined in Eq. (85) . with eigenvalue e k . This 
is equivalently stated by 



jw>(#5M$)- <B1) 



Decomposing this two-component equation into the top and 
bottom elements then reveals, directly, 

j rfr'L(r,r')Kit(r') + j dv'Miv, r>*(r') = e k u k (r), (B2) 

J" dr'M*(r,r')u k (r') + J" dr'L*(r,r')v k (r') = -e k v k (r). 

(B3) 

Taking the complex conjugates of the above equations then 
yields 

J" </r'L*(r,r>*(r') + J rfr'M*(r,r>*(r') = e * M *(r), (B4) 
j dr'M(r,r')u* k (r') + J rfr'L(r,r')v*(r')=-e t X(r). 

(B5) 

Combining Eq. (B4) and Eq. ( IB3b then yields a "left-hand" 
form of Eq. (BT) : 



dr(ut(r), -v*(r))X(r, r') = e;«(r'), -v*(r')). (B6) 



We now choose a normalization convention for the two- 
component eigenstates such that 



dr[\u k (r)\ 2 - Mr)! 2 ] = 1. 



(B7) 



Hence, applying Eq. ( IB6I ) onto a right eigenstate, where we 
can choose whether X(r, r') should act to the right [Eq. ( IB 1 1 11 
or the left [Eq. dB6H, reveals that 



JOT 



drdr'(u* k (r), -v*(r))£(r, r') ( ) = e =e t , ( B8 ) 



i.e., the eigenvalue is real. 

Thus, (HT(r), — vt(r)) is the corresponding left eigenstate, 
with eigenvalue e*, to the right eigenstate appearing in Eq. 

(EE). 

Furthermore, Eq. ( IB4b and Eq. ( IB3b imply that the complex 
conjugate of a right eigenstate is also a right eigenstate: 



/*'M$oH($'>) <B9) 



Now, in an equivalently manner to the derivation of Eq. 
from Eq. ( IB2b and Eq. dB5b we deduce that 



dr(-v k (r), u k (r))£(r, r') = -e k (-v k (r'), u k (r')), (BIO) 



i.e., that (— v k (r), u k (r)) is the corresponding left eigenstate, 
with eigenvalue -e k , to the right eigenstate appearing in Eq. 



As the eigenstates have different eigenvalues, they are or- 
thogonal, i.e., 



/ 
S 



dr[u* k ,(r)u k (r) - v|,(r)v*(r)] =6 kk ,, (Bl 1) 
dx\u v (r)Vi(r) - v k , (r)u k (r)] =0 (B 1 2) 



We note that setting u k (r) = 0(r) and v^(r) = on the one 
hand, and vt(r) = and wJ(r) = <p*(r) on the other, produces 
two eigenstates of eigenvalue zero. 

The identity can thus be decomposed as 



(B13) 



and, similarly, X(r, r') can be expressed as 
£(r,r0 ^("^) )("!(>•'), -vl(r')) 



(B14) 



This is the usual form of the spectral decompositions of the 
operator £(r, r') 145I1 . The two-component modes involving 
the condensate mode [which are also eigenstates of £(r, r')] 
do not explicitly appear as they have eigenvalue zero. 
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